This example is based on a set of local jpeg2000 images. Several parts have been commented out to minimize errors in the notebook. This component of the pipeline is the most modern. The previous version, on which many of the original reconstrctions were built, involved using ITK + c++ code to generate convert the points list into density stack. In theory this method should work on remote images from the API, but it has not been tested.
Sequence:
In [2]:
import sys, os
print 'Working in %s' % os.path.abspath(os.path.curdir)
# adding path to python modules
sys.path.append('../src/python')
We import the aibs module, a small wrapper created to query parts of the Allen Institute Data api api.brain-map.org easily
In [6]:
import aibs;
reload(aibs);
api = aibs.api()
Since this example uses local data, we manually define the specimen, raw file location, and other parameters.
Note There is a minor hack applied to this notebook to generate placeholder section images. The demo includes pre-transformed points to play with, but would break due to the missing raw images.
In [51]:
aSpecimen = {}
aSpecimen['subjectName'] = 'H08-0083.01.02'
#aSpecimen['location'] = '/vol/raw/UMB797/'
s = aibs.Specimen(initdata=aSpecimen, remote=True)
s.sectionRange = [71, 111] # we will only use sections from 19 to 129, inclusive
s.markersOfInterest = ['CALB1', 'RORB'] # no markers selected, use all of them
# uncomment as needed
#s.populateFromLocalImages() # searchs the location provided for images matching a known format
#s.getMarkerList(verbose=False) # compiles a list of possible markers to use
#all_images = s.getSectionImages() # filter all possible images by marrkers and section range
# minor hack to ONLY get this demo to work.
import glob
all_images = []
for rp in glob.glob('/data/reconstruction/specimens/H08-0083_01_02/register_points/*'):
a = aibs.SectionImage('-'.join(rp.split('-')[0:2]))
a.section_number = rp.split('-')[0]
a.metadata = {'height' : 0, 'width' : 0}
all_images.append(a)
We then import the processing module. This module was written to wrap various shell scripts & command line commands initially, but has since been expanded to include image processing steps itself, implemented in scikit image and numpy.
In [52]:
import pmip;
reload(pmip); # in case any development has occured since last import
# we create an instance of the class, passing it the experiment we defined above
pe = pmip.Processing(s)
# initializing the environment then creates the necessary directories for derived data to go
pe.initEnv();
# this is a utility command to see the total file counts in each directory
pe.listSubjectDirectory();
Import many of the modules we'll need downstream
In [53]:
import os, glob
import numpy as np
import scipy as sp
from scipy import ndimage, signal
import skimage
from skimage import color, filter, measure
from skimage.transform import pyramids
import workerpool
from mpl_toolkits.mplot3d import Axes3D
We then build several lists that will be used to interact with the points along the way.
In [54]:
#Collect the points list (centroids and areas), and transforms.
pe.collectImagesForGeneration()
generate_list = pe.processing_status['regpointc']
area_list = pe.processing_status['regpointa']
available_list = pe.processing_status['regcontrast']
xform_list = pe.processing_status['regxform']
points_to_make_with = pe.processing_status['regpoints']
for k in pe.processing_status.keys():
print k, len(pe.processing_status[k])
This transform was generated by first downsampling the native image by 2^4. The image was contrasted and padded to 2k x 2k from the center of the image. Then the image was shrunk by 50% to 1kx1k and registered using a centered 2d rigid registration
The original images were downsampled by a factor of 2 and run through a filter based point detection algorithm. This generated points in the downsampled space.
To get the transform into native space, take the center points of transform and rescale to native ( mult * 2^4)
e.g. 16 * cx, 16 * center.y
Next we need to transform the points detected in dsx1 space and place them into the 2kx2k*2^4 padded space. To do this, we will look up the section image dimensions from the specimen metadata and offset accordingly.
Thankfully, all of our images (of interest) are less than 32k*32k, meaning we only have to consider the possibility of insetting out points
In [55]:
def transformPoints(point_source, point_area, transform, imageInfo, lower_bound=30, upper_bound=100):
# print imageInfo.metadata
print imageInfo.metadata['path']
import numpy as np
from matplotlib.pyplot import scatter
from matplotlib import pyplot as plt
points = np.loadtxt(point_source, dtype=np.float64)
points_area = np.loadtxt(point_area, dtype=np.float64)
print point_area
#print points.shape
angle = cx = cy = dx = dy = 0
itk_fid = open(transform, 'rU')
read_pos = 0
for line in itk_fid:
items = line.split()
if read_pos == 3:
angle = float(items[1])
cx = float(items[2]) * 32
cy = float(items[3]) * 32
dx = float(items[4]) * 32
dy = float(items[5]) * 32
#print 'from this transform:'
#print '\tangle: %f' % angle
#print '\tcenter point: %f,%f' % (cx, cy)
#print '\tdelta: %f,%f' % (dx, dy)
read_pos +=1
left_offset = (32000 - imageInfo.metadata['width'])/2
top_offset = (32000 - imageInfo.metadata['height'])/2
new_points = []
olds = []
original_xform = []
for n,p in enumerate(points):
are = points_area[n]
if are >= lower_bound and are < upper_bound:
orx = (2*p[0])+ top_offset
ory = (2*p[1])+ left_offset
new_x = 2*p[0] + top_offset
new_y = 2*p[1] + left_offset
new_x -= cx
new_y -= cy
r_x = new_x*cos(angle) + new_y*-sin(angle);
r_y = new_x*sin(angle) + new_y*cos(angle);
good_x = r_x + cx - dy
good_y = (r_y + cy - dx);
newp = (good_x, good_y)
original_xform.append([orx, ory])
new_points.append(newp)
old_points = np.array(original_xform)
new_pointsss = np.array(new_points)
return new_pointsss
def getImageSizeForSectionNumber(number):
for a in all_images:
# print a.section_number
if a.section_number == number:
return a
This code shows how to transform the points using the above function and plot them via matplotlib
In [56]:
# #this will plot the first 2 pointlists
# import os
# def plotPointlist
# cnt = 0
# fig1 = plt.figure(figsize=(6,6), dpi=180)
# #ax = fig1.add_subplot(111, projection='3d')
# plt.xlim(0, 32000)
# plt.ylim(0, 32000)
# for ng, gen in enumerate(generate_list):
# if cnt < 2:
# basename = os.path.basename(gen).split('-')[0]
# for a in all_images:
# if a.section_number == int(basename):
# for n,avail in enumerate(available_list):
# si = os.path.basename(avail).split('-')[0]
# if basename == si:
# if n > 0:
# imageInfo = getImageSizeForSectionNumber(int(basename))
# points = transformPoints(gen, area_list[ng], xform_list[n-1], imageInfo, lower_bound=20, upper_bound=100)
# cval = cm.hot(cnt,1)
# zval = np.ones_like(points[:,0]) * cnt * .1
# if cnt == 0:
# y = plt.scatter(points[::2,1], points[::2,0], c='b', s=10);
# elif cnt == 1:
# y = plt.scatter(points[::2,1], points[::2,0], c='r', s=10);
# else:
# y = plt.scatter(points[:,1], points[:,0], c='g',s=10);
# cnt += 1
# plt.show()
Rather than apply the transformation each time we need registered points, we can do it once and save them to a new file. The code below does exactly that, but also restricts the output to points with an area between 20 to 100
In [71]:
# run the one-time conversion of points
def convertPoints():
import os
cnt = 0
for ng, gen in enumerate(generate_list):
basename = os.path.basename(gen).split('-')[0]
for a in all_images:
if a.section_number == int(basename):
for n,avail in enumerate(available_list):
si = os.path.basename(avail).split('-')[0]
if basename == si:
if n > 0:
imageInfo = getImageSizeForSectionNumber(int(basename))
points = transformPoints(gen, area_list[ng], xform_list[n-1], imageInfo, 20,100)
outputname = pe.dirs['regpoints'] + '/' + os.path.basename(gen) + '.reg'
import numpy as np
np.savetxt(outputname, points)
In [63]:
import numpy as np
rorb_count =0
calb1_count=0
for ppath in points_to_make_with:
# this corresponds to RORB labeling experiments
if '159713827' in ppath:
points = np.loadtxt(ppath, dtype=np.float64)
rescale = 16
softimg = np.zeros((32000/rescale,32000/rescale), dtype=np.float64)
for i,n in enumerate(points):
softimg[round(n[0]/rescale),round(n[1]/rescale)] += 1
improc = ndimage.filters.gaussian_filter(softimg, 200/rescale, mode='constant')
img_to_write = skimage.transform.pyramids.pyramid_reduce(improc, downscale=4)
ds4savename = os.path.join(pe.dirs['regdensity'], 'rorb-20-100-%04d.png' % (rorb_count))
sp.misc.imsave(ds4savename, img_to_write)
print ppath
print ds4savename
rorb_count+=1
# this corresponds to CALB1 labeling experiments
if '175707767' in ppath:
points = np.loadtxt(ppath, dtype=np.float64)
rescale = 16
softimg = np.zeros((32000/rescale,32000/rescale), dtype=np.float64)
for i,n in enumerate(points):
softimg[round(n[0]/rescale),round(n[1]/rescale)] += 1
improc = ndimage.filters.gaussian_filter(softimg, 200/rescale, mode='constant')
img_to_write = skimage.transform.pyramids.pyramid_reduce(improc, downscale=4)
ds4savename = os.path.join(pe.dirs['regdensity'], 'calb1-20-100-%04d.png' % (calb1_count))
sp.misc.imsave(ds4savename, img_to_write)
print ppath
print ds4savename
calb1_count+=1
Now that we have a registered set of points, we can generate our pseudo-density representation by convolving a 2D guassian against them. In other words, we're applying a gaussian filter to a binary image containing our points.
Once each density image is created, we can combine them into a single image stack and export it as one of many neuroimaging formats (shown here: nifti).
Note - This approach is not conservative and produces a volume that may appear correct, is a loose representation of what actual expression density is within the volume. Given that we're already producing a very sparse approximation of an entire volume from at most 10 serially-sectioned images, this approach is acceptable. When whole block imaging methods become available, this method will be obsolete.
In [75]:
# This will take each image from the RORB stack, create a numpy array, and export it as a NIFTI
import numpy as np
import nibabel as nib
rorb_count =0
calb1_count=0
volume_file = np.zeros((500,500,11), dtype=np.float64)
for ppath in points_to_make_with:
if '159713827' in ppath:
print ppath
points = np.loadtxt(ppath, dtype=np.float64)
rescale = 16
softimg = np.zeros((32000/rescale,32000/rescale), dtype=np.float64)
for i,n in enumerate(points):
softimg[round(n[0]/rescale),round(n[1]/rescale)] += 1
improc = ndimage.filters.gaussian_filter(softimg, 200/rescale, mode='constant')
img_to_write = skimage.transform.pyramids.pyramid_reduce(improc, downscale=4)
volume_file[:,:,rorb_count] = img_to_write
rorb_count+=1
img = nib.Nifti1Image(volume_file, np.eye(4))
stackname = os.path.join(pe.dirs['regstack'], 'rorb.nii.gz' )
img.to_filename(stackname)
In [15]:
# Same as above, except this time we will export it as an ANALYZE format, something easily readable by ImageJ
import numpy as np
import nibabel as nib
rorb_count =0
calb1_count=0
volume_file = np.zeros((500,500,15), dtype=np.float64)
for ppath in points_to_make_with:
if '175707767' in ppath:
points = np.loadtxt(ppath, dtype=np.float64)
rescale = 16
softimg = np.zeros((32000/rescale,32000/rescale), dtype=np.float64)
for i,n in enumerate(points):
softimg[round(n[0]/rescale),round(n[1]/rescale)] += 1
improc = ndimage.filters.gaussian_filter(softimg, 200/rescale, mode='constant')
img_to_write = skimage.transform.pyramids.pyramid_reduce(improc, downscale=4)
volume_file[:,:,calb1_count] = img_to_write
calb1_count+=1
#img = nib.Nifti1Image(volume_file, np.eye(4))
img = nib.AnalyzeImage(volume_file, np.eye(4))
stackname = os.path.join(pe.dirs['regstack'], 'calb1.img' )
img.to_filename(stackname)
We can plot one of our density images using matplotlib
In [79]:
fig = plt.figure(figsize=(16, 16), dpi=180)
ax = plt.imshow(improc, cmap=cm.gray)
cbar = fig.colorbar(ax)
imshow(improc, cmap=cm.gray)
plt.show()
In [74]:
# save a single density image
img_to_write = skimage.transform.pyramids.pyramid_reduce(improc, downscale=4)
ds4savename = os.path.join(pe.dirs['regdensity'], os.path.basename(ppath) + '.png')
sp.misc.imsave(ds4savename, img_to_write)
print ds4savename
Finally, if you made it this far, we can also explore the points list in a bit more detail
In [65]:
# for curiousity, show histrogram of point size distribution
pointdist = '/data/reconstruction/specimens/H08-0083_01_02/detect_points/029-175707767-DSx1.jpg.area'
points = np.loadtxt(pointdist, dtype=np.float64);
def reject_outliers(data, m=2):
return data[abs(data - np.mean(data)) < m * np.std(data)]
points = reject_outliers(points)
plt.hist(points, bins=8000)
plt.xlim(0, 400)
Out[65]:
In [66]:
# taking a look at the histogram bins & values
histoval = np.histogram(points, bins=100)
histoval
Out[66]:
In [67]:
# let's fit a curve to it!
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a*np.exp(-b*x) + c
x_vals = histoval[1][0:100]
y_vals = histoval[0][0:100]
popt, pcov = curve_fit(func, x_vals, y_vals)
print popt
print x_vals.max()
In [68]:
# and plot our curve
x = np.linspace(1,x_vals.max(), 1000)
_y = np.zeros_like(x)
for n,_x in enumerate(x):
_y[n] = func(_x, popt[0], popt[1], popt[2])
plt.scatter(x, _y)
plt.xlim(1,250)
Out[68]:
In [69]:
# Subtracting curve
fix_y = np.zeros_like(y_vals)
for n,nx in enumerate(x_vals):
fix_y[n] = y_vals[n] - func(nx, popt[0], popt[1], popt[2])
plt.scatter(x_vals, fix_y)
plt.xlim(0,250)
Out[69]:
In [ ]: